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ABSTRACT 

We present a fully Lagrangian conservation form of the general relativistic 
hydrodynamic equations for perfect fluids with artificial viscosity in a given arbi- 
trary background spacetime. This conservation formulation is achieved by choos- 
ing suitable Lagrangian time evolution variables, from which the generic fluid 
variables rest-mass density, 3-velocity, and thermodynamic pressure have to be 
determined. We present the corresponding equations for an ideal gas and show 
the existence and uniqueness of the solution. On the basis of the Lagrangian 
formulation we have developed a three-dimensional general relativistic Smoothed 
Particle Hydrodynamics (SPH) code using the standard SPH formalism as known 
from non-relativistic fluid dynamics. One-dimensional simulations of a shock tube 
and a wall shock are presented. With our method we can model ultra-relativistic 
fluid flows including shocks with relativistic 7-factors of even 1000. 



Subject headings: hydrodynamics — methods: numerical — relativity — shock 
waves 



INTRODUCTION 



Modeling ultra-relativistic fluid flows is a great challenge for any relativistic hydro code. 
Numerical difficulties arise from strong relativistic shocks and from narrow physical structures 
(Norman &: Winkler 1986). Typical examples in astrophysics for such extreme conditions are 
proto-stellar jets and blast waves of supernovae explosions. In recent years, the development 
of numerical algorithms for relativistic fluid dynamics went mainly along two different lines. 
First, there are the so called High Resolution Shock Capturing (HRSC) methods, which allow 
to obtain numerically discontinuous solutions of the relativistic hydrodynamic equations by 
solving local Riemann problems between adjacent numerical cells. Some recently developed 
HRSC techniques are those of Font et al. (1994), Falle & Komissarov (1996), Romero et 
al. (1996), Banyuls et al. (1997), Wen, Panaitescu, & Laguna (1997), Pons et al. (1998), 
and Komissarov (1999). However, employing analytic solutions of the Riemann shock tube 
problem, it is not surprising that these HRSC codes produce almost exact numerical solu- 
tions of flow structures including discontinuities. In the second type of algorithms, shocks are 
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treated numerically not as fluid discontinuities, but are rather spread over some small length 
with the help of an artificial viscosity. These algorithms either solve the dynamical equations 
of relativistic hydrodynamics on an Eulerian grid such as the finite difference schemes of 
Hawley, Smarr, &; Wilson (1984b) and Norman &; Winkler (1986) or by using computational 
nodes following the fluid motion such as the Lagrangian Smoothed Particle Hydrodynamics 
(SPH) methods of Kheyfets, Miller, & Zurek (1990) and Laguna, Miller, k Zurek (1993). All 
the latter methods seem to be limited even for mildly relativistic flows containing shocks. 
Norman & Winkler (1986) suggest that the appearance of numerical inaccuracies and in- 
stabilities is due to the time derivative of the relativistic 7-factor in the energy equation of 
these numerical schemes. This additional time derivative of a hydrodynamic variable renders 
the system of evolution equations non-conscrvativc. For non-relativistic hydrodynamics, it 
is well-known that a numerical method based on non-conservative equations can produce a 
solution which looks reasonable but is entirely wrong if shocks or other discontinuities are 
involved (see LeVeque 1997 for a corresponding analysis of Burger's equation). 

The main intention of this paper is to present a conservation formulation of the relativistic 
hydrodynamic equations that is designed for numerical methods. This particular formulation 
allows hydro codes to resolve ultra-relativistic shocks numerically with relativistic 7-factors 
of even 1000 by means of an artificial viscosity rather than using a Riemann solver. The 
conservation form of the relativistic equations is obtained by choosing suitable Lagrangian 
variables. Unfortunately, these Lagrangian dynamical variables can be expressed in terms 
of the generic fluid variables rest-mass density, 3-velocity, and thermodynamic pressure only 
through a set of nonlinear algebraic equations. However, we show that these equations can 
be solved analytically in a unique way. 

For numerical work we discretize our set of partial differential equations by the Smoothed 
Particle Hydrodynamics (SPH) method introduced by Lucy (1977) and Gingold &; Monaghan 
(1977). In recent years, SPH became a popular computational tool for numerically modeling 
complex three-dimensional fluid flows in astrophysics. This is primarily due to its computa- 
tional simplicity and the absence of a computational grid. Furthermore, SPH is adaptive in 
the sense that its computational nodes follow the fluid. 

SPH has been already applied to relativistic fluid flows. Kheyfets et al. (1990) developed 
a relativistically covariant version of the SPH technique by modeling the contact interac- 
tions with spatial smoothing functions in the local comoving frame of the fluid. Laguna et 
al. (1993) applied the SPH method to the so called ADM formalism of general relativity due 
to Arnowitt, Deser, & Misner (1962). In their relativistic SPH formulation they modifled the 
flat space kernels of the Newtonian SPH method which can become anisotropic and are then 
no longer invariant under translations which leads to additional terms in the SPH equations. 
In contrast to these methods, we have developed a fully three-dimensional general relativistic 
SPH code on the basis of our set of Lagrangian conservation equations. This code employs 
the standard SPH approach as used in Newtonian theory with spherically symmetric ker- 
nels for all particles. Our code is restricted to ideal fluids with artificial viscosity using the 
ideal-gas equation of state. The influence of the fluid on the spacetime metric is neglected. 
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therefore, we consider only a background spacetime with a given but otherwise arbitrary 
metric. In addition, we neglect the fluid's self-gravity and do not account for radiation or 
electromagnetic effects. 

The layout of this paper is as follows. In section ^ we derive the formulation of the rela- 
tivistic hydrodynamic equations in Lagrangian conservative form and present the equations 
for calculating the generic fluid properties rest-mass density, 3-velocity, and thermodynamic 
pressure from our suitably chosen Lagrangian variables. Section ^ gives a review of the stan- 
dard SPH method that we use, including a prescription of its application to relativistic fluid 
flows in curved spacetimes and the implementation of an artificial viscosity. Numerical test 
calculations of one-dimensional ultra-relativistic flows with discontinuities are presented in 
section 0. Throughout this paper we set the speed of light c and the Boltzmann's constant 
k to unity, i.e., c = k = 1. Latin indices ■ ■ ■} run from 1 to 3, Greek indices {//, i', ■ ■ ■} 
from to 3; the signature of the metric is (— , +, +, +). 



2. THE LAGRANGIAN CONSERVATION FORMULATION OF THE 
GENERAL RELATIVISTIC IDEAL FLUID EQUATIONS 

Our starting point is the covariant formulation of the equations of relativistic hydrody- 
namics for a perfect fluid with artificial viscosity: the local conservation of baryon number 

ipu''),, = (1) 
and the local conservation of energy-momentum 

T^^^ = (2) 

with the stress-energy tensor of a perfect fluid 

Tf"^ = {pw + q)u''u'' + ip + q)g^'' . (3) 

The semicolon denotes the covariant derivative and the Einstein summation convention is 
used. Here, p is the rest-mass density, the 4-velocity of the fluid, p the isotropic ther- 
modynamic pressure, w = l+ e + p/p the relativistic specific enthalpy with specific internal 
energy e, q the artificial viscous pressure, and g^'^ the spacetime metric. All thermodynamic 
quantities in the stress-energy tensor ^ are measured in the local rest frame of the fluid. 
The artificial viscous pressure q, which appears in equation can be introduced into the 
stress-energy tensor either by the substitution p ^ p + q or, equivalently, from the stress- 
energy tensor of a viscous fiuid (Misner, Thorne, &: Wheeler 1973; Landau & Lifshitz 1991) 
by ignoring the shear viscosity and heat conduction and replacing the bulk viscous pressure 
by q. An explicit expression for q in terms of velocity gradients will be given in the subsequent 
section. 

For a Lagrangian formulation of relativistic hydrodynamics suitable for SPH one has 
to break the unity of time and space inherent in the covariant formulation. This can be 
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accomplished by applying the ADM formalism of Arnowitt et al. (1962), where the spacetime 
is decomposed into an infinite foliation of spatial hypersurfaces of constant coordinate 
time t by writing the line element as 

ds'^ = g^^dx^'dx" = - (a^ _ pi^.^j ^^2 ^ 2/3idx' + r]ijdx'dx^ . 

Here, a is the lapse function, /?* the shift vector, and rjij the spatial metric induced on Et with 
Pi = 'HijP'' 'Hii'n^'' = ^i'' ■ III this paper we consider only given background spacetimes, 
i.e., we do not solve the Einstein equations. Thus, g^^, ct, f3^, and r]ij are given analytic 
functions of both space and time. From the definitions of a and /3* the basis vector field dt of 
the coordinate basis {dt,di} can be decomposed into normal and parallel components with 
respect to the hypersurfaces Sj 

dt = an + P'di , (4) 

where n is the unit time-like vector field normal to the slices Sj, i.e., n • di = 0. Observers 
having n as 4-velocity are at rest in the slices St — they are called Eulerian observers. In 
the basis {n, 5j} the 4-velocity of a fiuid has the representation 

whereas in the coordinate basis {9t,9j} it follows from equation 

with 

= av' - . 

From the normalization condition of the 4-velocity u'^u^ = — 1 the relativistic 7- factor is 
given by 

1 

7 



^1 - rjijV^yj 

In the following, we use both, the 3-velocity of the fluid v'^ measured by Eulerian observers 
and the 3-velocity in the coordinate basis. 

We now derive the Lagrangian equations of relativistic hydrodynamics where the La- 
grangian or total time derivative d/dt is defined as 

d (y. 

— = -uf'd^ = dt + v'di . 
dt ^ 

From the law of baryon-number conservation (||) we obtain 

= i^^D-9y, (5) 
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where g is the determinant of the spacetime metric g^'^ , and the relativistic rest-mass density 
D* is defined by 



D* 



-9-P = V^IP 



a 



(6) 



Here, rj is the determinant of the spatial metric r]ij obeying y/rj = ^J—gja. With the above 
definition of D* the relativistic continuity equation (P) has the same form as in the non- 
relativistic case. Note, however, that in equation (|5|) the expression for diV^ is not the diver- 
gence of a 3- vector on but rather a sum of partial derivatives. Eulerian observers measure 
the relativistic rest-mass density T) = —pu • n = p^. By rewriting equation (^) for D, we 
obtain the relativistic continuity equation 

which, in contrast to equation (|5|), contains an additional source term. In this paper we will 
use both definitions of a relativistic rest-mass density, i.e., D and D* = ^/rjD. 

In order to derive our relativistic expressions for the energy and momentum equations, 
we first re-express the conservation law of energy-momentum (|2|) by using the continuity 
equation (||) 







pu 



w -\ — I n 
P 



(7) 



One can easily show that the covariant derivative in the first term of equation (|^ can be 
written as 



u 



w -\ — \ u 



7 d 
a dt 



w -\ — ] u 



2\ p) ^"^'''^ 



Thus, equation becomes 

w + - \ Ui, 



d_ 

di 



P 



a 

PI 
1 

15* 



d^i{p + q)- - {pw + q)u°u'^gaf3,i 



dfi [V^{P + q)] 



(8) 



where we have used the identity 



--d„ 



~9) 



1 



9 ^9a/3,^, 



Taking the spatial components of equation (^), i.e., p = i, and using 

Ui = gi^uf" = (3i - + rjij-v^ = jTjijv^ , 
a a 



we get the momentum equation 
d 



dt 



5,; 



1 

D* 



di [V^{p + q)] 



9al3 



(9) 
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with the relativistic specific momentum Si defined by 



(10) 



The component DSi = —T{n,di) of the stress-energy tensor field T is the relativistic mo- 
mentum density in the z-direction measured by Eulerian observers. The momentum equation 
(|9|), which, in a similar form, was also used by Laguna et al. (1993), can be applied immedi- 
ately to the SPH method because it contains only spatial derivatives on its right hand side. 
In the Eulerian formulation an expression similar to equation (^) without artificial viscosity 
is given by Hawley, Smarr, & Wilson (1984a) for the momentum variable DSi, which is more 
convenient for a Eulerian description. The Newtonian limit yielding the non-relativistic Euler 
equation is obvious from equation (P). 

Taking the /i = 0-component of equation (^), we obtain the relativistic energy equation 



d 
di 



w + - ] uo 
P 



1 

15* 



dt 



-g{p + q)] 



-g 



(11) 



which, unfortunately, has time derivatives of hydrodynamic variables on both sides. Re- 
expressing the left hand side of equation (|TT1) as 



w + - \ Uq 
P. 



P 



with 



2\ 7 , ^ 7 i 
I - + Pi-v' 
a a 



a 



and rewriting the first term on the right hand side of equation (O) using equation IM 



D 



-dt [^/^{p + q)\ 



Id 

[V^{P + Q)] - Jy^di + q)] 



d_ 

di 



-g- 



-p + q 
D* 



1 

'd* 



di 



-g{p + q)v' 



we obtain the relativistic energy equation 



d_ 
lit 



aE - (3'Si 



1 

IT* 



di 



-gip + qW 



+ 



hgap 



with the total relativistic specific energy E defined as 

p + q 



E 



li; + - I 7 ■ 



D 



(12) 



(13) 



The component DE = T{n,n) of the stress-energy tensor field T is the total relativistic 
energy density measured by Eulerian observers. As in the case of the relativistic momentum 



equation (|9|) the right hand side of equation (12) contains no time derivatives of hydrodynamic 
variables. It is, therefore, well suited for the SPH method in contrast to the energy equation 
used by Hawley et al. (1984a, 1984b) and Laguna et al. (1993) which has a non-conservative 
form containing two total time derivatives of hydrodynamical variables separately on both 
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sides of the equation. Norman & Winkler (1986) suggest that the additional time derivative 
of the relativistic 7-factor in the energy equation gives rise to numerical inaccuracies and 
instabilities. Note that for the special relativistic case equation (|l2|) without an artificial 
viscosity can also be found in Monaghan (1992). In the non-relativistic case, equation ( p!^ ) 
yields the energy equation for the total non-relativistic specific energy \v\^ /2 + em written in 



a form similar to our expression (12) 



dt 



\v\ + [ WN + 



2 V PN J PN 



■^3. 

PN 



where the index N denotes Newtonian quantities and wn = £n +p/ Pn is the non-relativistic 
specific enthalpy. 

To close our system of hydrodynamical equations (|5|) , (|9|) , and (^) , we have to add an 
equation of state of the form p = p(p, e), which relates the thermodynamic pressure p to the 
rest-mass density p and the specific internal energy e. We restrict ourselves to the ideal-gas 
equation of state given by 

p={T- l)pe (14) 

with the ideal-gas adiabatic constant T. 

Our system of ideal fluid equations is now complete. We have derived the relativistic 
hydrodynamic equations (^), (^, and (12) in Lagrangian conservative form similar to their 



Newtonian counterparts. This was achieved by choosing suitable hydrodynamic variables D*, 
Si , and E defined in equations (^) , ( [To| ) , and ( |T^ ) . Because of the equation of state ( p^ ) and 
the use of the 3- velocity for moving particles in the Lagrangian numerical methods, we now 
need to calculate the generic hydrodynamic quantities p, u*, and p from these variables by 
solving a highly nonlinear system of equations. If an artificial viscous pressure q is included, 
a severe difficulty arises from q being usually expressed in terms of velocity gradients. Since 
there is no time evolution equation for q, the character of the dynamic equations is thus 
changed which is the usual situation in the hydrodynamics of viscous flows. Due to the cou- 
pling of and q, the suitably chosen dynamic variables actually have to be solved iteratively 
for the generic hydrodynamic variables. However, artificial viscosity operates only in the 
vicinity of shock transitions and is zero everywhere else. For simplicity one can calculate 
the generic variables in a single iteration taking q from the previous time step. Using the 
expression w = 1 +p/{Gp) with G = 1 — l/F for the total relativistic specific enthalpy, the 
specific energy E from equation (|l3|) can be written as 



G / 1 ~ 

E = — iw — \) h 7 — 

7 V UP 

--\w + - + (-f . (15) 

7/ 7 V 7/ P 

Solving equation ( [T5| ) for the relativistic specific enthalpy w and adding the term q/p, we 
obtain 

q Ej-G 
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where the variable E is given by 

E = E + ^ . 
TD 

Using 

S^ = rj'^SiS,= (w + ^y{j'-l) (17) 

and inserting expression (16), the relativistic 7- factor can be determined explicitly from 

= - E^^ 7^ + 2GEj^ + (e"^ - 2GS^ - G^) 7^ - 2GE-f + {l + S^) (18) 

as the root of a polynomial of degree four. In the appendix we show that a solution of 
equation (|l^) exists and that it is unique for all allowed values of G, E, and S"^. With the 
value of 7 known, one can calculate first the rest-mass density p from equation (|6|), then the 
thermodynamic pressure p from equation (^) and the equation of state ( |l^ , and finally the 
velocity v"^ from equation (10) using r]^^ Sj = {w + q/p)'yv'^. 



3. THE SPH EQUATIONS 

In this section we derive the SPH formalism for our set of relativistic hydrodynamic 
equations (|5|) , (^) , ( |l^ ) , and ([l^) . Since in the previous section we have obtained the dynamic 
equations in an appropriate Lagrangian form, we can proceed in a way that is completely 
analogous to the Newtonian case. Before we derive our relativistic SPH equations, we give 
a brief introduction into the standard SPH formalism, which was invented independently 
by Lucy (1977) and Gingold & Monaghan (1977). For a review of the SPH method see 
Monaghan (1992). 

SPH is a numerical method which discretizes the dynamic equations on a set of nodes, 
called particles, moving with the fluid. The final discrete equations are obtained in two 
separate steps. First, all hydrodynamic functions on T,t are smoothed over a certain volume 
with the help of a smoothing kernel Ty(|r — r'|, h), i.e., for a continuous function f{r) one 
has 

f{r) = J f{r') W{\r - r'\,h) dr' + 0{h'^) , (19) 

where r = {x*} is the set of spatial coordinates x*. The kernel is a smooth (differentiable) 
function with compact support of size /i, the so called smoothing length, it is normalized to 
unity 

j W{\r -r'\,h)dr' = 1 , 

and consequently, the smoothing error in equation ( ]l9| ) is 0{h?). The above integrals extend 
over a region on around r that contains the support of W . An example for W is the 
cubic-spline kernel of Monaghan & Lattanzio (1985), and in our simulations we use this 
kernel throughout. The second step consists of an approximate evaluation of the above 
integral (pjt) at the particle positions 

I f{r')W{\ra-r'\,h)dr' ^Y.-^-b=--U)a ■ (20) 
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Here, a and b label the particles which are distributed in space with number density n(r). 
We define fa = f{ra) for any function / and Wab = W{\ra — r-fel, h). The important point is 
that, applying the smoothing and discretization operations ([l9|), (pO|), it is possible to derive 
approximate expressions for derivatives 

=E-^-^afe ' (9^f) =T.-fb-'^-Wab , (21) 
^ rib \ I a ^ Ub 

where / = {/*}, V = {di}, and '^aWab is the gradient of the kernel — taken 

with respect to the coordinates of particle a. 

With the discretization method described above we are now in the position to derive the 
SPH equations for our system of relativistic hydrodynamic equations (|5|) , @ , (|l^) , and (|lj) , 
and nothing has to be changed compared to the standard (non-relativistic) SPH method. In 
the equations (|2^) and (^T|) we replace the number density Ua by the mass per particle ma, 
i.e., Ua = D^/ma- Dropping the angle brackets ( ) from now on, equation (pO|) leads to the 
SPH representation for the relativistic rest-mass density 

Dl = Y^mbWab . (22) 
b 

From equation (^ and using 

D*V ■ V = V-{D*v) - V • VD* , 

we obtain the SPH form of the relativistic continuity equation 

j^D: = - ^ rUbiVb - VaYVaWab • (23) 
b 

Here, v = {v'} is just a shorthand notation for the set of components v^, whereas V • v 
stands for diV^ and is not the divergence of a vector field v. As a consequence, the total mass 
is conserved exactly. Applying the Lagrangian time derivative to equation (22) with the 



smoothing length being constant in space and time, one can show that the SPH expression of 



the relativistic rest-mass density D* in equation (22) automatically satisfies the relativistic 
continuity equation. Thus, in SPH the relativistic rest-mass density D* can be computed 
either from equation or equation (|^). If, however, the density D = jp is used instead of 
D*, a modification of the standard SPH method is required. Therefore, Laguna et al. (1993) 
multiply the flat-space kernel W{\r — r'\,h) of Newtonian SPH with l/y/r}{r'). In a curved 
spacetime this factor makes the kernel anisotropic, violates its translation invariance, and 
leads to additional terms in the SPH approximation of derivatives. Applying the relation 
^* — \/^-^' relativistic SPH formulation of Laguna et al. (1993) can be cast into the 
standard SPH scheme which is considerably simpler. We therefore suggest that there is no 
need to modify SPH for given background spacetimes if the continuity equation is used in 
the conservative form (||) without source terms. 
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in order to derive the SPH form of the relativistic momentum equation we start from 
equation (P). Rewriting the pressure gradient as 



-9 



P + Q 
D* 



D 



*2 



P + 9„ 
+ 

D* 



-a 



we obtain 



dt 



Sa 



, (Pa^ 



Pa + qa Pb + Qb 
2 + n*2 



D* 



^aWab 



(Pa + qa)Va (In - -T^^V a (^a/?). 



(24) 



where Sa = and the metric gradients V In y/—g and ^Qa/s can be calculated analyti- 

cally. Thus, only the pressure gradient term in equation ( p^ has to be smoothed, and it has 
been symmetrized to conserve linear and angular momentum in SPH. Next, we proceed with 
the relativistic energy equation ([T^). In the expression 



1 



D* 



we replace the first term on the right hand side by 



l-V-[{p + q)v] = l-v.V{p + q) + ^V-v 



p + q\ , P + q 



+ 
+ 



1 



D* 



2 

p + q 

D 



D* 



V • V 



P + Q 
D* 



*2 



[V-{D*v) - V • VD* 



With this combination of terms we obtain the SPH representation of the relativistic energy 
equation 



It 



aE - 13' S, 



^/-ga fPa + Qa , Pb + 



D 



*2 



{Va + Vb)-VaWab 



D* 



(Pa + qa)Va-Va (in + -T^^{9aP,t)a 



(25) 



Note that equation ( [25| ) and the momentum equation (^4|) contain identical symmetric factors 
in front of the kernel gradients. As outlined in section |2|, this form of the relativistic energy 
equation is well suited for SPH because it contains no time derivatives of hydrodynamic 
variables on its right hand side. Again, there is no need to smooth the derivatives of the 
metric 'V\n^/^^ and gaf3,t- The last equation that needs to be considered is the equation of 
state (14), and its SPH formulation reads 



Pa = (r - l)pa£a 



(26) 



- 11 - 



This completes our set of relativistic SPH equations (22), (24), (25), and (26) for computing 
general relativistic fluid flows in given arbitrary background spacetimes. 

We now focus our attention to the implementation of an artificial viscosity term which 
is necessary to handle shock fronts. For the use of an artificial viscosity in the standard SPH 
method see Monaghan & Gingold (1983). In our simulations, we have used the following 
artificial viscous pressure 




-aCahaiV ■ v)^ + phliV ■ v)l if (V • ^), < 



otherwise 



with 



(V 



\rab? + ehl^ 



(27) 



(28) 



where Cq = yjTpa/ {paWa) is the relativistic sound velocity measured in the rest frame of the 
fluid, d, /? and e are numerical parameters, Vab, ^ab stand for the differences Vab = Va — I'fe, 
fab = fa — fb, and hab = {ha + hb)/2 is the mean value of the smoothing lengths of particles a 
and b. Without the enthalpy Wa, the expression ( p7| ) for q is almost equivalent to the standard 
SPH form of the artificial viscosity invented by Monaghan & Gingold (1983). Including Wa 
into equation (27), the parameters a and (3 can be chosen to be of order unity even for shocks 
with ultra-relativistic values of 7. The d-term, which is linear in the velocity differences, is 
similar to a physical shear and bulk viscosity, and the quadratic /3-term is the standard von 
Neumann-Richtmyer (1950) artificial viscosity used in finite difference methods for handling 
high Mach-number shocks. According to Monaghan & Gingold (1983), the un-smoothed 
representation of the velocity divergence in equation ( p8|) acts more directly on the relative 
motion of particle pairs and leads to a damping of irregular oscillations in shock transitions. 
In equation (|^), the parameter e has been introduced to avoid singularities, and a typical 
value is e = 0.1. 



Since the expression (|27| ) for the artificial viscosity is not Lorentz invariant, we also 
performed calculations with a relativistically covariant formulation of the artificial viscous 
pressure 



Pa 





where 



i^'';.)a 



1 



an 



-aCahaOa + 



d-fa 

dt 



2 

a" a 



if 0a < 

otherwise 
^)a + ;|(ln^/^)„ 



(29) 



However, this expression contains a time derivative of the 7-factor which destroys the explicit 
nature of the Lagrangian form of the hydrodynamic equations. One possibility to circumvent 
this problem is to take the backward time difference approximation [7a (t) — "fa{t — At)] /At 



with time step At for the time derivative dja/dt. Since the non-covariant form (27) turned 
out to be quite appropriate for the resolution of shock structures, we used this approach in 



all our simulations and did not pursue the covariant relation (29) 
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To improve the local resolution of SPH, we allow the smoothing length h to vary in space 
according to the relativistic rest-mass density D* via 



K = {ho 



ma 



D* 



l/d 



(30) 



where d denotes the dimension of the spatial slices T,f For particles of equal mass the 
scaling law (|30|) indicates that the number of particles within the support of the kernel W 
is approximately kept constant in time. Thus, in regions where the gas is compressed, the 
smoothing length is increased while it is decreased in rarefaction zones. Since the smoothing 
length is now a function of the density D*, equation (|22|) (or eq. |^3|) is now a nonlinear 
implicit relation for the density. We solve this approximately by inserting D* from the 
previous time step into equation (^) to obtain an estimate for ha- Since the smoothing 
length is now a function of position and time, the SPH form of the continuity equation 
contains additional terms (Monaghan 1992). However, no such terms appear if the relativistic 
rest-mass density D* is calculated from the computationally simpler equation (p^). 



4. NUMERICAL TESTS 

Although we have developed a fully three-dimensional general relativistic SPH code, we 
restrict ourselves in this paper to the standard analytic test bed of one-dimensional special 
relativistic shock problems: the shock tube and the wall shock. For each simulation we 
show four diagrams of the numerical results together with the analytic solution for the fluid 
variables 3- velocity v := = v^, rest-mass density p, thermodynamic pressure p, and specific 
internal energy e. These variables are functions of the coordinate x := G [0,100]. The 
quality of the simulations is measured in terms of the relative error with respect to the 
analytic solution, i.e., for each hydrodynamic function / an error Af is calculated from 

= TTTo — '^Ifb- fb\ ) 

-'■^ J max f^—i 

where the sum is over all particles and /" stands for the analytic solution which has the 
maximum value /^ax- 



a) Shock Tube Tests 

First we consider the shock tube problem, where initially a fluid at rest is divided by a 
diaphragm into two regions of different densities and internal energies. When the diaphragm 
is removed, a rarefaction wave travels into the warm and dense medium and a compression 
wave into the cold and lower density fluid. Between the two media a so called contact 
discontinuity is present. For the analytic solution of the relativistic shock tube problem we 
refer to Taub (1948), McKee & Colgate (1973), Hawley et al. (1984a), and Marti & Miiller 
(1994). 



- 13 - 



3 10 



2 10 



1 10 



20 40 60 80 100 
x-coordinate 




40 60 80 
x-coordinate 



100 



0.75 



2 0.5 



0.25 




3 10 



2.5 10 



2 10 



20 40 60 80 100 
x-coordinate 



1.5 10 



20 40 60 80 100 
x-coordinate 



Fig. 1. — Numerical result of an SPH calculation with 1000 particles (open circles) and 
analytic solution (solid line) for the non-relativistic shock tube problem (note that the units 
are arbitrary except for the velocity which is measured in units of c). The intermediate 
pressure, velocity, and relativistic 7-factor are pm = 0.303, Vm = 2.93 x 10~^, and 7^ = 
1.000004, the positions of the shock, the contact discontinuity, and the head and tail of the 
rarefaction wave are Xg = 83.2, Xc = 67.6, = 27.6, and xt = 48.7, respectively, and the 
velocity and relativistic 7-factor of the shock are Vg = 5.54 X 10-^ and 7^ = 1.000015. 



As any relativistic hydro code has to be tested for the Newtonian limit, we start our series 
of shock tube simulations with a non-relativistic test case. Figure || shows the corresponding 
numerical and analytic solution for a gas with F = 1.4 at time t = 6000 with the initial 
conditions p = 10^, e = 2.5 x 10"^ for x < 50, and p = 0.125 x 10^, e = 2 x 10~^ for x > 50 
(note that the units are arbitrary except for the velocity which is measured in units of c) . The 
numerical calculation was performed with 1000 particles, initial smoothing length h = 0.6 
(« 10 interactions per particle), and artificial viscosity parameters a = 1 and [5 = 2. In the 
initial distribution, particles were placed on a uniform grid, and the particles at x > 50 have 
a mass ten times smaller than the particles at x < 50. In the calculation of Figure [l| the 
largest relative error occurs in the fluid velocity where Aw = 1.1%. 
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Fig. 2. — Numerical result of an SPH calculation with 1000 particles (open circles) and ana- 
lytic solution (solid line) for the relativistic shock tube problem. The intermediate pressure, 
velocity, and relativistic 7-factor are pm = 1-45, Vm = 0.714, and 7m = 1-4, the positions 
of the shock, the contact discontinuity, and the head and tail of the rarefaction wave are 
Xs = 87.3, Xc = 82.1, Xh = 17.8, and xt = 57.5, respectively, and the velocity and relativistic 
7-factor of the shock are Vg = 0.828 and 7s = 1.8. 



Next, a mildly relativistic shock tube is investigated with initial conditions p = 10, 
e = 2 {x < 50), and p = 1, e = {x > 50). Figure |2| shows the numerical result and 

the analytic solution for a gas with F = 5/3 at time t = 45. We used the same numerical 
parameters as in the non-relativistic shock tube problem, i.e., 1000 particles, h = 0.6, a = 1, 
and (3 = 2. The largest relative error is Av = 1.0%. As can be seen in Figure ||, the velocity 
profile of a relativistic rarefaction wave is no longer linear as in the Newtonian case because 
of the relativistic velocity addition formula. Comparing our results with the simulations of 
Hawley et al. (1984b) and Laguna et al. (1993) for the same 7 = 1.4 shock tube, we note 
that numerical inaccuracies due to the non-conservative formulation of their energy equation 
clearly show up in their results for the relativistic rest-mass density and specific internal 
energy. To investigate the convergence properties of our numerical method, we performed a 
calculation of the 7 = 1.4 shock tube with 10000 particles and initial smoothing length h = 0.1 
(~ 20 interactions per particle). As can be seen in Figure ^, the numerical calculation covers 
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Fig. 3. — Numerical result of an SPH calculation with 10000 particles (points) and analytic 
solution (solid line) for the relativistic shock tube problem. The intermediate pressure, ve- 
locity, and relativistic 7-factor arc p,n = 1.45, = 0.714, and 7„j = 1.4, the positions of the 
shock, the contact discontinuity, and the head and tail of the rarefaction wave are Xg = 87.3, 
Xc = 82.1, Xh = 17.8, and xt = 57.5, respectively, and the velocity and relativistic 7-factor of 
the shock are Vg = 0.828 and 7^ = 1.8. 



the analytic solution almost exactly, and the largest relative error is reduced to Av = 0.2%. 

When the relativistic 7-factor of the shock tube is increased by increasing the initial 

specific internal energy ratio, the region between the leading shock front and the trailing 
contact discontinuity becomes extremely thin and dense. Thus, without specially designed 
adaptive methods, it will be impossible to resolve these Lorentz contracted shells of matter, 
which are typical for relativistic fluid flows. In addition, the speciflc internal energy at the 
contact discontinuity may become negative in the SPH simulations. In order to avoid these 
difficulties, we consider now a simplified problem of a single shock front without the presence 
of a contact discontinuity and a nonlinear rarefaction wave. 
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Fig. 4. — Numerical result of an SPH calculation with about 250 particles (open circles) 
and analytic solution (solid line) for the relativistic 7j = 1.8 wall shock problem. The post 
shock properties are pp = 5.66, pp = 1.52, and Ep = 0.803, and the position and velocity of 
the shock front are Xg = 64.3 and Vs = —0.178 with 7^ = 1.02. 



h) Wall Shock Tests 

As a second test of our numerical method we modeled the wall shock problem of a cold 
relativistically moving fluid flowing towards a solid wall. As the fluid hits the wall, a shock 
front forms, which then travels upstream against the incoming fluid producing a hot and 
dense post shock region of zero velocity. 

Figure ^ shows the numerical result and the analytic solution of a mildly relativistic wall 
shock for a gas with F = 4/3 at time t = 200 moving to the right with the reflecting wall at 
X = 100. The uniform initial fluid properties are Di = 1, Vi = 0.832 with 7^ = 1.8, and Si = 
10^^. Initially, particles of equal mass are uniformly distributed in the simulation domain. 
At the location of the solid wall we use reflecting boundary conditions. The simulation of 
Figure |^ showing about 250 particles was performed with the initial smoothing length h = 3 
(~ 10 interactions per particle) and artificial viscosity parameters a = 0.25 and (3 = 0.5. The 
rest-mass density p, the thermodynamic pressure p, and the specific internal energy e show 
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Fig. 5. — Numerical result of an SPH calculation with about 250 particles (open circles) and 
analytic solution (solid line) for the ultra-relativistic 7^ = 1000 wall shock problem. The post 
shock properties are pp f« 4, pp 4/3 x 10^, and Sp ^ 1000, and the position and velocity of 
the shock front are Xg ~ 33.3 and Vg ~ —1/3 with 7^ ~ 1.06. 



a spike-like feature (see Fig. ^), which is known in the literature as "wall heating" (Norman 
& Winkler 1986). The largest relative error appears in the specific internal energy where 
Ae = 1.1%. 

The results of an ultra-relativistic shock simulation are shown in Figure ^. The initial 
velocity Vi is increased to Vi = 0.9999995 which corresponds to a relativistic 7-factor of 
7j = 1000. All other parameters are identical with those of the previous wall shock calculation. 
Neglecting the pre-shock specific internal energy ej , one can show that in the ultra-relativistic 
limit — > 1 the post shock properties pp, pp, Ep, and the shock velocity Vg are given by 
Pp = D,jT /{T — 1), Pp = TDiji, Ep = 7j, and Vg = — (F — 1). For the wall shock problem of 
Figure |5| we thus obtain pp A, pp 4/3 x 10'^, Ep 1000, and Vg —1/3. At the time 
t = 200 the shock front has moved the distance \vg\t « 66.7 to the left reaching the spatial 
position Xg « 33.3. The largest relative error in this case is Ap = 1.1%. 

A calculation of a mildly relativistic wall shock problem (with 7j = 2.24) using artificial 
viscosity was also performed by Hawley et al. (1984b) with their time-explicit Eulerian finite 
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Fig. 6. — Numerical result of a finite difference calculation with a grid of 250 zones (open 
circles) and analytic solution (solid line) for the ultra-relativistic 7^ = 1000 wall shock prob- 
lem. The post shock properties are pp ^ A, pp ^ 4/3 x 10^, and Sp ^ 1000, and the position 
and velocity of the shock front are Xs ~ 33.3 and Vs ~ —1/3 with 7^ 1.06. 



difference code. They tried several formulations of their non-conservative energy equation by 
omitting and including various artificial viscous pressure terms. They found that the qdtj 
term in their energy equation is unstable even for mildly relativistic wall shocks. This is 
verified by the calculations of Norman & Winkler (1986) using an implicit adaptive-mesh 
finite difference method, which shows the requirement of an implicit technique to handle the 
qdt'j term. The fact that our method is able to treat even ultra-relativistic shocks without 
problems is entirely based on our conservation formulation of the relativistic hydrodynamic 
equations where no additional time derivatives of hydrodynamic variables appear such as the 
qdt'j term. In order to verify that this numerical stability is not restricted to the SPH method, 
we also performed numerical simulations of the wall shock problem with a simple explicit finite 
difference upwind scheme. Figure ^ shows the numerical result and the analytic solution for 
the 7j = 1000 wall shock at time t = 200 using a grid of 250 zones, and the maximum relative 
error occurs for the specific internal energy where Ae = 2.2%. 
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5. SUMMARY 

We have derived a fully Lagrangian conservative form of the general relativistic equations 
of hydrodynamics for a perfect fluid with artificial viscosity in a given arbitrary background 
spacetime. This has been achieved by choosing suitable Lagrangian time evolution variables. 
These variables are connected to the generic fluid variables rest-mass density p, 3- velocity v*, 
and thermodynamic pressure p through a set of nonlinear algebraic equations. For an ideal 
gas, we have shown that these equations can be reduced to a single fourth-order equation 
with a unique solution which can be explicitly calculated in terms of various roots. For more 
complex equations of state the solution of the above algebraic equations has to be performed 
numerically, and the question of uniqueness is more complicated. Using our Lagrangian 
formulation, we have developed a three-dimensional general relativistic SPH code based on 
the standard SPH approach of non-relativistic fluid dynamics. The important point is that 
all metric factors from covariant derivatives have been absorbed in the definition of the 
Lagrangian variables and thus no longer appear in the fluid equations. As a result, the SPH 
kernels remain spherically symmetric and are the same for all particles. This is an essential 
difference to the covariant SPH approach of Kheyfets et al. (1990) using kernels defined in the 
local comoving frame of the fluid, and to the relativistic SPH method of Laguna et al. (1993) 
where, in curved spacetime, the kernels are anisotropic and have no translational symmetry 
which leads to additional terms in the SPH equations. The relativistic continuity equation 
(HI), for example, contains no source term, hence we can identify the relativistic rest-mass 
density D* = y/rj^p as the appropriate SPH density for smoothing the equations. 

In this paper we have restricted ourselves to the numerical simulation of two different 
one-dimensional examples, i.e., the special relativistic shock tube and the wall shock. An 
empirical error estimate is obtained from a comparison of the numerical results with the 
corresponding analytic solutions. The SPH calculations with 1000 particles show a typical 
maximum relative error of about 1%, and an increase of the particle number with a decreasing 
smoothing length reduces this error in a uniform way. The wall-shock problem can be solved 
without any numerical difficulties for very large 7- factors (at least 7 = 1000). The shock-tube 
case suffers from a resolution problem if 7 is large because a thin shell of matter builds up 
between the shock front and the contact discontinuity. This can only be resolved with the 
help of additional adaptive methods which increase the number of particles in this zone. 

An important ingredient in our SPH formulation is the treatment of shock structures by 
an artificial viscosity rather than using a Riemann solver. This is very easy to implement 
even for the two- or three-dimensional case. Introducing an artificial viscous pressure q is 
considered as a purely numerical tool that acts as a filter to smear out steep gradients in the 
hydrodynamic functions and suppresses unphysical oscillations. Thus, there is no need to use 
a covariant expression for q in terms of velocity derivatives as long as the time evolution and 
the jump conditions of shocks are represented correctly. Since the jump conditions follow 
from the conservation of mass, momentum, and energy across the shock, it is obvious that 
a conservative form of the relativistic hydrodynamic equations is important for handling 
discontinuities in relativistic fiuid dynamics numerically. The simulation of relativistic fiows 
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with large 7-factors including shocks is the traditional domain of HRSC methods. However, 
with our formulation of the general relativistic hydrodynamic equations it is possible to model 
such flows using an artificial viscosity independent of the underlying numerical method, i.e., 
SPH or finite difference schemes. 
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A. APPENDIX 

In section |2| we have derived equation (18) from which the relativistic 7- factor can be 
calculated analytically as the root of a polynomial of degree four for given artificial viscous 



pressures q. Equation ( [18| ) is restricted to an ideal-gas equation of state (14). We will now 
show the existence and uniqueness of the solution of equation (^) for all allowed values of 
G, and E. 

First, we rewrite equation (^) by substituting 7 = 1 + 5, 5 G [0, cxd) and obtain 







E"^] 5^ + 2[2S' - 2E' + GE]5- 



a4 



as 



+ (gS^ - 5E^ + 6GE - 2GS^ - G^^ 6^ 
V . ■' 

012 

+ 2 (2S^ -E^ + 2GE - 2GS^ - G^) 5 + 5^(1 - G) 



(Al) 



ao 



With the degree of freedom / > 3 the ideal-gas adiabatic constant T = 1 + 2/ f and the 
variable G = 1 — l/F lie in the range 



i< r < 



< G < 



3 ' - 5 

From equations (16) and ( [l7| ) the variables E and 5^ are given by 



(A2) 



E 



G 



7 



G 



7 ] w -\ 



c2 ~2 
D = W 



7 



(7' - 1) > 



(A3) 



where we have defined w = w + q/ p > 1. Using the relations (A2) and the expressions (^ 
we obtain for the coefficient 04 in equation (Al) 



04 



S^-E^ 



1 -2G + 



G^ 
7^ 



w 



2 - 2G 1 



G_ 



G2 



w 



< . 



(A4) 



Thus, 5^ is limited to the range 



0< <E^ 



a) Existence 



With the coefficients 04 < and ag 
positive root of 5. 



S /V > equation (Al) has at minimum one 
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b) Uniqueness 



To show the uniqueness of the solutions of equation (^) for 5, we investigate the changes 
of signs of the coefficients as, 02, and ai depending on the variable 5^: 

1 



03 < ^ <E'^ 



02 < ^ < 



-GE =: S- 



1 



ai < <^ < 



2(3 - G) 
1 

2(1 -G) 



3 ! 



5E^ - 6GE + G' 



S? 



2 ' 



and 



E-G 



-.Sf 



Using 



E-G = {w 
and the relations ([A2|), one can show that 

Sf > , 

Si - sf 



7 



+ 7- G > 



r 



s'. 



> 



< 



i + 2r 

TG^ 

i + 2r 
1 

2(3- G) 
1 

2(3- G) 
5 



2 - r)^2 _^ 2(r - 1)G^ - TG^ 
[(2 - r) + 2(r - 1) - r] = , and 



6GE + G^ 



2(3 - G) 

The relations (^, (^, and (^ lead to 



(^2 _ < E^ - GE < si . 



6GE + GS 



(A5) 



(A6) 



(A7) 



o<sf< si < si 



(A8) 



i) For ao = S^/F^ > or S"^ > 0, respectively, we obtain from the relations ( |A4| ) and 
(|A8|) the following table of signs for the coefficients a^, 03, 02, ai, and ao depending on S^: 



S^ 


04 03 


a2 


ai 


ao 


si <S^<E^ 


- + 


+ 


+ 


+ 


52 = Si 


- 


+ 


+ 


+ 


Si <S^< Si 




+ 


+ 


+ 


52 = Si 







+ 


+ 


sj <s^< si 






+ 


+ 


52 = 52 









+ 


< 52 < ^2 








+ 



For all values of 52 G (0,^2) ^j^e series of coefficients 04, 03, 02, ai, and ao has exactly one 
change of sign. Therefore, in 6 £ (0, 00) there exists only one root of equation (Al) for 6. 
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ii) For ao = S"^ /T^ 



the relations (A4) and 
(5 = is the only positive root of equation (| 



or = 0, respectively, equation ( |Al| ) has the root 5 = 0. With 
the coefficients 04, 03, a2, and oi are all less than zero. Thus, 
for 5. 



To summarize, with the restriction of the ideal-gas equation of state ( [T^ ) a solution of 
equation (Al) for 5 or of equation (|T^) for 7, respectively, exists and is unique. Thus, the fluid 
variables rest-mass density 3- velocity -u*, and thermodynamic pressure p can be calculated 
analytically in a unique way from the variables D* , Si, E, and q from equations (|T8[), (^), 
(0), (0), and (0). 
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